function [corrEEG] = eegArtifactCorrection(eegData,eogData,intervalIdx)

% =========================================================================
% [corrEEG] = eegArtifactCorrection(eegData,eogData,intervalIdx)
%           corrects the data by the amount of
%           eye-motion recorded by the VEOG and HEOG. 
%           eegData, eogData - the timeseries data nCh X nTime X nTRials
%           intervalIdx - time samples for estimation
% =========================================================================
% 
% Purpose: activity due to eye-motion has to be removed
%
% INPUT:   3-dimensional dataset with the format    
%          channel-by-timepoints-by-trialnumber
%
% OUTPUT:  a new matrix with the same format as the input but corrected by
%          the estimated amount of activity due to eye-motion
%
% REQUIREMENTS: The statistics toolbox

% 100301 SD wrote it
% 100309 Revision: Aeye will be estimated by the solution of the
%           regressions
%           SD wrote it
% 100312 Revision: VEOG is channel 28
%           SD wrote it
%---------------------------------------------------
% eye-motion removal according to PARRA et al (2005)

if nargin<3,
    intervalIdx=true(1,size(eegData,2));
end

corrEEG = zeros(size(eegData));

    for k = 1:size( eegData,3 );
        crit        = eegData(:,intervalIdx,k)';
        
        %-------------------------------------------------
        % Set of regression equations will be set up
        
        equatSet    = [ones(length(intervalIdx),1) eogData(:,intervalIdx,k)'];
        
        %-------------------------------------------------
        % yields a 29-by-3 matrix of regression parameters
        
        Aeye        = (equatSet\crit)';
        
   
        % reconstruction of eeg signal not affected by artifacts (8)
        I               = eye( size( Aeye,1 ) ); 
        corrEEG(:,:,k)  = ( I - Aeye * pinv( Aeye ) ) * eegData(:,:,k);
    end
end